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Abstract 

The  power  spectral  density  of  frequency  fluctuations  of  an  oscillator  is  gen¬ 
erally  modeled  as  a  sum  of  power  laws  with  integer  exponents  (from  -2  to  -f-2). 
However,  a  power  law  with  a  fractional  exponent  may  exist.  We  propose  a  method 
for  measuring  such  a  noise  level  and  determining  the  probability  density  of  the 
exponent.  This  yields  a  criterion  for  compatibility  with  an  integer  exponent. 
This  method  is  based  upon  a  Bayesian  approach  called  the  reference  analysis  of 
Bernardo- Berger.  The  application  to  a  sequence  of  frequency  measurement  from 
a  quartz  oscillator  illustrates  this  paper. 


INTRODUCTION 

It  is  commonly  assumed  that  Sy{f),  the  power  spectral  density  (PSD)  of  frequency 
deviation  of  an  oscillator,  may  be  modeled  as  the  sum  of  5  power  laws,  defining  5 
types  of  noise: 

+2 

Syif)  =  Y  (1) 

a=-2 

where  is  the  level  of  the  /"  noise.  But  it  may  be  noticed  that  models  with  non-integer 
exponents  are  occasionally  used. 

The  estimation  of  the  noise  levels  is  mainly  achieved  by  using  the  Allan  variance  [i], 
which  is  defined  versus  the  integration  time  t  as; 

(2) 
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In  the  frequency  domain,  the  Allan  variance  may  be  considered  as  a  filter.  If  the  Allan 
variance  versus  the  integration  time  r  is  plotted,  the  graph  exhibits  different  slopes, 
each  slope  corresponding  to  a  type  of  noise; 

(r^(r)  =  <1^  S'y(/)  =  and  q  = -p  -  1.  (3) 

The  estimation  of  yields  an  estimation  of  the  noise  level  ha. 

However,  this  curve  may  exhibit  an  exponent  ^  which  seems  to  be  non-integer.  Does 
this  mean  that  the  corresponding  PSD  is  not  compatible  with  the  5  power  law  model? 
In  this  paper,  we  propose  a  method  for  estimating  the  most  probable  value  of  this 
exponent  in  order  to  solve  this  ambiguity.  This  method  is  applied  to  an  example  of 
stability  measurement. 


CLASSICAL  STABILITY  ANALYSIS  OF  AN  OSCILLATOR 

Sequence  of  frequency  measures 

Figure  1  shows  average  frequency  measures  F/,  of  a  10  MHz  quartz  oscillator  compared 
to  a  cesium  clock.  The  sampling  rate  is  10  s  and  the  integration  time  of  each  frequency 
measure  is  also  10  s  (sampling  without  dead  time). 

In  order  to  obtain  dimensionless  y),  samples,  we  must  subtract  the  nominal  frequency 
i/Q  (10  MHz)  from  the  frequency  measures  and  normalize  by  i/qI 


Variance  analysis 

Figure  2  is  a  log-log  plot  of  the  Allan  deviation  of  the  quartz  samples  versus  the 
integration  time  r.  A  least  squares  fit  of  these  variance  measures  (solid  line),  weighted 
by  their  uncertainties,  detects  only  two  types  of  noise:  a  white  noise  and  an  noise. 
The  corresponding  noise  level  estimations  are: 

ho  =  (2.2  ±  0.4)  ■  10~^s  at  1<t  (68%  confidence) 

=  (2-3  ±  0.6)  •  at  Icr  (68%  confidence) 

(for  the  assessment  of  the  ha  noise  levels  and  their  uncertainties,  we  used  the  multi¬ 
variance  method  described  in  M). 

However,  for  large  t  values  (corresponding  to  low  frequencies),  the  variance  measures 
move  away  from  the  fitted  curve.  Two  explanations  are  possible: 

•  instead  of  an  f~^  noise,  there  is  a  noise  whose  non-integer  exponent  is  contained 
between  -2  and  -3  ; 

•  since  the  uncertainty  domains  of  the  variance  measures  contain  the  fitted  curve, 
this  apparent  divergence  may  be  due  to  a  statistical  effect. 


In  order  to  choose  between  these  two  explanations,  we  decided  to  estimate  the  proba¬ 
bility  density  of  the  exponent  with  a  Bayesian  approach. 


BAYESIAN  APPROACH 

Principle 

The  goal  of  all  measurement  is  the  estimation  of  an  unknown  quantity  9  from  measures 
i.e.  determining  p{6\(,)^  the  density  of  probability  of  the  quantity  B  knowing  the 
measures  The  Bayesian  theory  is  based  on  the  following  equality  W: 

P(<9|f)  ocp(^|6i)7r(6>)  (5) 

where  is  the  distribution  of  the  measures  ^  for  a  fixed  value  of  the  quantity  9  and 

Tr(9)  is  the  a  priori  density  of  probability  of  the  quantity  B,  i.e.  before  performing  any 
measurement. 

The  determination  of  this  a  priori  density,  called  the  prior,  is  generally  one  of  the  main 
difficulties  of  this  approach  (particularly  in  the  case  of  total  lack  of  knowledge!  ).  In 
this  paper,  we  use  the  .leffrey’s  prior  which  ensures  properties  such  as  invariance  M. 


Spectral  density  and  covariance  matrix 

Let  us  define  the  vector  y  whose  components  are  the  N  samples.  We  assume  that  y 
is  a  Gaussian  vector.  The  probability  distribution  of  y  is: 


exp 

(271)^  vie 


(6) 


where  C  is  the  covariance  matrix.  Since  Sy{f)  is  the  Fourier  transform  of  Tlyir),  the 
autocorrelation  function  of  the  frequency  deviation,  the  general  term  of  C  is: 


=  2  y  Sy{f)  cos  (27r/(ti  -  tj))  (if. 


(7) 


Equation  (7)  reveals  the  key  role  played  by  the  spectral  density  of  the  noise  in  the 
expected  fluctuation.  We  will  present  a  general  method  for  estimating  the  parameters 
of  the  model  for  Sy{f). 


Assumed  model  for  the  spectral  density 

We  assume  that  the  sequence  of  frequency  measures  is  composed  of  a  white  noise  yw, 
whose  variance  (i.e.  the  level)  is  unitary,  and  of  a  red  noise  whose  level  is  unknown, 
multiplied  by  the  real  standard  deviation  of  the  white  noise  cTyj  (cr„  is  easily  estimated 
from  high  sampling  rate  frequency  measurement): 


y  =  {yw  +  yr)<rw  ■  (8) 

This  yields  the  following  model  for  Sy(f): 

Syif)  =  ho  +  ha  ■  f°‘  with  —  3  <  q  <  —  1  (9) 

where  ho  —  2.2  •  10“®s,  ha  and  a  are  the  unknown  parameters. 

Let  us  denote  y„,  the  normalized  vector: 


J/n  =  S/w  +  Vr- 


(10) 
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The  corresponding  normalized  PSD  S’„(/)  is; 

Sr.if)^l  +  H r  (11) 

where  is  an  amplitude  factor  whose  meaning  will  be  explained  below  (see  equation 
(23)). 

Statistical  model 

The  part  of  the  spectral  density  due  to  the  red  noise  yr  may  be  written: 

Srif)  =  H-Uc-r.  (12) 

We  used  the  Bernardo- Berger  analysis  D;  4]  for  estimating  the  unknown  parameter 

e= 

-  Construction  of  the  estimators: 


Let  us  introduce  the  orthonormal  basis  of  {po, . .  .,pj, . . .  defined  such  as  the 

i^^  component  of  pj  is: 

Pij  =  (13) 

where  ti  is  the  date  of  the  frequency  measure  and  pj{t)  is  a  polynomial  of  degree  j, 
satisfying  the  orthonormality  condition  l^l: 

JV-l 

Y^PjiU)  ■Pk{ti)  =  ^jk-  (14) 

i=0 

It  can  be  shown  than  the  scalar  product  of  a  vector  pj  by  the  noise  vector  y  is  an 
estimate  of  the  noise  spectrum  for  a  given  frequency  fj  Let  us  denote  such  an 
estimate  applied  to  the  normalized  noise: 

i]-Pj-yn-  (15) 

Practically,  we  limited  to  16  the  number  of  estimators  pj  (from  degrees  0  to  15)  for  lim¬ 
iting  the  computation  and  because  the  high  degrees,  estimating  the  high  frequencies, 
are  less  informative  for  a  red  noise. 

Moreover,  in  order  to  ensure  convergence  for  very  low  frequencies  (even  if  the  low  cut¬ 
off  frequency  tends  towards  0),  the  polynomials  must  satisfy  the  moment  condition  [5>  6]; 
the  minimum  degree  jmin  of  a  polynomial  to  ensure  convergence  up  to  an  exponent  a 
is:  -a-1 

jmin  >  - ■  (16) 

Since  we  have  assumed  a  <  —3,  the  first  2  estimators  (po  and  pi)  must  be  removed. 
Thus  we  have  rj  =  14  estimators  {p2, . .  .,pi5}  and  n  =  14  estimates  {6,  •  •  •  ,4i5}- 

-  Construction  of  the  priors; 

The  covariance  matrix  defined  in  relationship  (7)  is  an  ensemble  average  of  the  different 
■estimate  products  over  an  infinite  number  of  realizations  of  this  process; 

C  =  (17) 

^ij  ~  '  ^j)  ■  (18) 
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As  for  the  noise  vector  yn,  the  estimate  vector  ^  may  be  split  into  two  terms,  according 
to  equations  (10)  and  (15): 

i  =  (19) 

The  covariance  matrix  may  also  be  split: 

c  =  +  {irC)  =  ln+H-u^-  Via)  (20) 

where  /„  is  the  identity  matrix  in  3?"  and  the  general  term  of  the  matrix  T(a)  is: 

[b(a)L  =  2  /  /“  cos  {2TTf{U  -  ij))  df.  (21) 

Jl/T 

The  high  cut-off  frequency  fh  in  (21)  is  the  Nyquist  frequency  and  T  is  the  total 
duration  of  the  sequence. 

Let  e,(n)  denote  the  eigenvector  of  V{a)  and  7,(0)  its  eigenvalue  (i  €  -  1}). 

The  averaged  quadratic  norm  of  the  estimate  vector  ^  is: 

n  — 1 

=  n  +  F  ■  ^7.(«)  =  (ll^^ll')  +  (ll^rll')  ■  (22) 

iz:.D 


The  factor  u 
and 


,  is  chosen  in  such  a  way  that,  for  F  =  1,  the  averaged  quadratic  norms 
11611^^  are  equal: 

(23) 


En—  I 

i=Q 


(u) 


The  direct  problem  is  now  solved  since  f  is  a  vector  of  3?”  with  a  probability  distribution 
given  the  parameter  6  equal  to: 

p(^i^)  =  ~ - L^exp(-ie‘C-'0-  (24) 

^  (27r)"/2v1q 

Denoting  “Tr(M)”  the  trace  of  a  matrix  M  and  X  the  matrix  defined  as: 

X  =  Ua  ■  V(a)  (25) 


the  Fisher  information  matrix  /(&)  is  (see  W): 


(C-ifC-ig) 
FTr(C“iXC-i^^ 


da  > 


FTr(C-iXC-i£) 
Tr  (C-^XC“iX) 


(26) 


The  Jeffrey’s  prior  7r(d)  is  defined  as: 

-(6)  =  vW)J-  (27) 

The  parameter  ^  is  a  two-dimensional  parameter  composed  of  the  exponent  parameter 
a  and  of  the  amplitude  parameter  F.  Since  we  are  mostly  interested  in  a,  H  is  called 
a  nuisance  parameter. 

In  presence  of  nuisance  parameter,  Bernardo  and  Berger  suggested  that  a  should  first 
be  fixed  and  the  conditional  prior  7:{H\a)  computed  for  that  value.  The  full  prior  is 
then: 

7i-(0)  =  7v{H\a)  ■  7r(a)-  (28) 
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The  conditional  prior  -K{H\a)  is  given  by: 


=  VWm^\  (29) 

where  [I{0)]i2  =  [■^(^)]2i  and  [/(^)]22  are  the  elements  of  the  Fisher  information 

matrix  J{0). 

The  prior  for  a  may  be  computed  as: 

7r(a)  =  c  •  exp  7r(if  |q)  In  H)f'^  d>^  (30) 

where  c  is  a  normalization  coefficient  ensuring  that  /  7r(a)da  =  1  and: 

=  (31) 

This  prior  is  plotted  in  Figure  3. 

-  Construction  of  the  posteriors: 


According  to  the  Bayes  theorem,  the  posterior  probability  distribution  is  given  by: 

p(^|^)7r(6») 


’  Jp{m7r{o')de'- 

The  posterior  probability  distribution  for  a  is  then  given  by: 

f  p(4  la,  If  )7:(ffla)7:(a)dE 


P(a|0  = 


/  /p(^1q;',  If')7r(If'lQ')7r(a')dIf'da' ' 


(32) 


(33) 


RESULTS  AND  DISCUSSION 

Compatibility  with  an  integer  exponent 

Figure  4  shows  the  posterior  probability  distribution  for  the  exponent  a  of  the  red 
noise  using  the  Bernardo-Berger  prior. 

The  exponent  value  obtained  for  the  maximum  of  likelihood,  just  as  for  the  maximum 
of  the  distribution,  is  a  =  -2.2. 

However,  a  =  -2  is  fully  compatible  with  this  prior  distribution.  Thus,  we  may  conclude 
that  the  apparent  divergence  between  the  variance  measures  and  the  fitted  curve  in 
Figure  2  is  probably  due  to  a  statistical  bias  of  the  data.  The  spectral  density  Syif)  is 
then  compatible  with  the  following  model: 

Syif)=ho+h.2r--  (34) 

Noise  level  estimation 

Selecting  an  exponent  value  a  —  -2,  we  obtained  the  posterior  probability  distribution 
plotted  in  Figure  5.  As  in  the  variance  analysis,  we  chose  a  confidence  interval  of  68% 
(16%  probability  that  /i_2  is  smaller  than  the  low  bound  and  16%  probability  that  ft_2 
is  greater  than  the  high  bound): 

h-2  =  ^2.3  ^  ^  at  la  (68%  confidence) 
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The  difference  between  the  maximum  likelihood  value  (/i_2  —  2.2991  ■  10  and  the 

variance  analysis  value  (/j_2  —  2.2949  •  10~^“s’'' )  is  only  0.18%. 

However,  the  confidence  intervals  given  by  these  two  methods  are  quiet  different.  The 
main  difference  concerns  the  symmetry  of  the  variance  analysis  interval;  in  this  case, 
we  don’t  take  into  account  the  fact  that  the  noise  levels  are  positive,  whereas  the  prior 
of  the  Bayesian  approach  is  null  for  negative  values  of  /i_2. 

Moreover,  the  variance  analysis  interval  seems  to  be  a  bit  underestimated. 

CONCLUSION 

The  variance  analysis  is  an  useful  tool  for  a  quick  estimation  of  the  noise  levels  in 
the  output  signal  of  an  oscillator.  However,  a  negative  estimate  of  a  noise  level  may 
occur.  Generally,  in  this  case,  this  value  is  rejected  and  the  corresponding  noise  level 
is  assumed  to  be  null.  On  the  other  hand,  although  the  Bayesian  method  is  a  bit 
heavier,  it  takes  into  account  properly  the  a  priori  information,  and  gives  a  more 
reliable  estimation  of  these  noise  levels  and  especially  of  their  confidence  intervals. 

However,  the  main  advantage  of  the  Bayesian  method  concerns  the  verification  of  the 
validity  of  the  power  law  model  of  spectral  density.  Each  time  the  model  is  suspected, 
such  an  approach  should  be  used  in  order  to  estimate  the  exponent  of  the  power  law. 
In  particular,  this  method  should  be  very  interesting  for  the  study  of  the  and 
noise,  whose  origin  remains  mysterious  M. 
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Figure  1 :  Sequence  of  frequency  measures 
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Figure  2:  Allan  deviation  of  the  sequence  of  frequency  measures 


Figure  4:  Posterior  probability  density  for  the  power  q 


